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Abstract. - Dynamics of many-body Hamiltonian systems with long range interactions is studied, 
in the context of the so called a— HMF model. Building on the analogy with the related mean field 
model, we construct stationary states of the a— HMF model for which the spatial organization 
satisfies a fractional equation. At variance, the microscopic dynamics turns out to be regular and 
explicitly known. As a consequence, dynamical regularity is achieved at the price of strong spatial 
complexity, namely a microscopic inhomogeneity which locally displays scale invariance. 
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' The out-of-equilibrium dynamics of many body systems 
subject to long range couplings defines a fascinating field 
of investigations, which can potentially impact different 
domains of applications [T] . In a long range system every 
constituent is simultaneously solicited by the ensemble of 
microscopic actors, resulting in a complex dynamical pic- 
ture. This latter scenario applies to a vast realm of fun- 
damental problems, including gravity and plasma physics, 
and also extends to a rich variety of cross disciplinary stud- 
ies [3]. In particular, and with specific emphasis on the pe- 
culiar non equilibrium features, long range systems often 
display a slow relaxation to equilibrium. They are in fact 
trapped in long-lasting out of equilibrium regimes, termed 
ip the literature Quasi Stationary States (QSS) which bear 
distinct characteristics, as compared to the corresponding 
deputed equilibrium configuration. A paradigmatic repre- 
sentative of long range interactions, sharing the mean field 
viewpoint, is the so called Hamiltonian Mean Field model, 
which describes the coupled evolution of TV rotators, pop- 
ulating the unitary circle and interacting via a cosines like 
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potential. In the limit of infinite system size the discrete 
HMF model is described by the Vlasov equation for the 
evolution of the single particle distribution function [5]. 
This leads to a statistically based treatment, inspired by 
the seminal work of Lynden-Bell, which revealed the ex- 
istence of two different classes of QSS, spatially homoge- 
neous or inhomogeneous [4l[5|. More recently, and still 
with reference to the HMF case study, stationary states 
have been constructed using a dynamical scheme which 
exploited the formal analogy with a set of uncoupled pen- 
dula [S]. This represented a substantial leap forward in 
the understanding of the dynamical properties of the QSS 
in the HMF model, beyond the aforementioned statistical 
approach. Indeed, it was understood that the microscopic 
dynamics in the inhomogeneous stationary state is regular 
and explicitly known, an observation which contributed to 
shed light onto the puzzling abundance of emerging regu- 
lar orbits as revealed in [7]. 

These last results have been though obtained in the 
framework of mean field systems: The actual distance 
between particles does no explicitly appear in the HMF 
potential. In this letter we aim at bridging this gap, by 
focusing on the so called a— HMF model [5]. We shall 
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here concentrate on the long range version of the model 
which implies assuming < a < 1 In this long range 
version the model behaves at equilibrium as the HMF, 
see for instance [Tl[9l[Tni. We may then ask if the same 
correspondence applies to the out of equilibrium dynam- 
ics. QSS exists for the a— HMF model, as shown in [TT] . 
However, can one still appreciate the asymptotic trend to- 
wards regularity? How does the spatial organization im- 
pact on the aforementioned features? It was also shown 
recently that fractional calculus may be a crucial ingre- 
dient when dealing with long range systems |12| . and we 
shall see how this point enters the picture in the consid- 
ered case. To anticipate our findings, we shall here show 
that all stationary states of the HMF model are shared by 
the a— HMF model: Particles still exhibit regular orbits, 
at the price of a enhanced microscopic spatial complexity, 
which materializes in a small scale inhomogeneity being 
locally scale invariant. Let us start by introducing the 
governing Hamiltonian which can be cast in the form: 



TV 



pi 

2 



1 

2N 



^ 1 - cos(g, - qj) 



(1) 



where qi stands for the orientation of the rotor occupying 
the lattice position i, while pi labels the conjugate momen- 
tum. The quantity ||i — j|| denotes the shortest distance 
on the circle of circumference N. The coupling constant 
between classical rotators decays as a power law of the 
sites distance. The HMF limit is recovered for a = 0. For 
N even, we have 
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which guarantees extensivity of the system. The equa- 
tions of motion of element i are derived from the above 
Hamiltonian and can be written as follows 
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where use has been made of the following global quantities 
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These identify the two components of a non-local mag- 
netization per site, with modulus Mi = y^Cf + Sf , and 
phase (fii — arctan(5'i/Ci). In doing so, one brings into ev- 
idence the formal analogy with the HMF setting. Indeed 
we notice that each individual a— HMF particle obeys a 
dynamical equation which closely resembles that of a pen- 
dulum. This observation represented the starting point of 
the analysis carried out in [5] , where stationary states were 
constructed from first principles. More specifically, the au- 
thors of [B] imagined that the system of coupled rotators 



reached a given stationary state, characterized by a con- 
stant magnetization in the limit for TV — >■ oo . Then, the 
HMF model is mapped into a set of N uncoupled pendula, 
constrained to collectively return a global magnetization 
identical to that driving their individual dynamics. Self- 
consistency is hence a crucial ingredient explicitly accom- 
modated for the formulation proposed in [S]. Technically, 
the stationary state is built by exploiting the ergodic mea- 
sure on the torus originating from the pendulum motion, 
a working ansatz that we cannot invoke in the context 
of the a— HMF, due to site localization. Eventually we 
will overcome this difficulty by considering the continuous 
limit. We notice that, for large N , and since < a < 1, 
the following estimate applies 
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We can then use expression ([SJ in Eq.(|3]) and, as cxo 
, introduce the continuous variables x = i/N and y = j /N 
and arrive to the following Riemann integral 
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where ||a; — y|| represents the minimal distance on a circle 
of circumference one. By invoking the Riemann-Liouville 
formalism on the circle, we recognize the fractional inte- 
gral /i~". To be more precise, the Riemann-Liouville for- 
malism on a circle cannot be defined in a consistent way, 
as the Riemann-Liouville operations cannot map periodic 
functions into periodic functions. However, the integra- 
tion is performed over the whole circle in Eq.Q. And, 
in fact, we are relying on fractional integrals of the po- 
tential type, which means we are actually considering a 
linear combination of Weyl fractional derivatives, see for 
instance jl3j for more details. Consequently we can write 

C(^) = ^r(l-a)/i-"(cosQ(.x)) . (8) 

A similar relation holds for the S{x) component. Studying 
the a— HMF dynamics implies characterizing the evolution 
of the scalar fields q{x,t) and p{x,t) which are ruled by 
the fractional (non-local) partial differential equations 



dt 
dp 



p{x,t) 

^r{p) (-sin(g)/^ (cosq) -^cos(g)/^ (sing)) . 



where /x = 1 — a. 

At variance with the simpler HMF {a ~ 0) model, the 
spatial organization q{x) matters within the general set- 
ting < a < 1, an observation which, as anticipated 
above, poses technical problems to a straightforward ex- 
tension of the analysis in [5] . 

In order to carry on the study and mimick the mean 
field situation, we turn to considering the particular con- 
dition where C'{x) = (C) = constant, where (• • •) denotes 
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a spatial average. At the same time, and without losing 
generality, we require /'^(sing) = 0, which corresponds 
to setting to zero the S component of the magnetization. 
Hence, we choose to deal with a constant magnetization 
amount M such that C{x) = M. 

Note that the translation invariance along the lattice 
is also likely to statistically lead to the configuration 
here scrutinized |10] . For finite size systems, the iden- 
tity C{x) — M can be solely matched by the trivial state 
where all q arc equal. Conversely, in the infinite N limit, 
the latter assumption, or cquivalcntly dC/dx — 0, implies: 
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where we recalled expression ([5]) and where the opera- 
tor stands for the fractional derivative. Trivial states 
(q(x)=Cte) as evidenced in the finite size approximation 
are also solutions of this equation. However, as we will 
argue in the following, the continuum limit enables us to 
compute an independent set of admissible solutions. This 
task is accomplished by exploiting the fact that a < 1, 
which makes the integral J l/||a;||"fix convergent near 
and that the function l/\\x — y\\" is smooth. Noticing 
that we rewrite Eq. ^ as: 

using the regularity ofl/||a; — we can extract it from 
the integral: 
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with yk G [k/L,(k + l)/L] The above approximation can 
be rigorously derived via a detailed expansion that we en- 
close in the annexed appendix. Taking advantage from the 
latter expression and aiming at reproducing a non trivial 
state with C{x) = M, we now assign the q values on a cir- 
cle so that in any small interval the average of cosg(x) is 
constant and equal to the global, constant magnetization 
M of the system. This procedure implies a peculiar spa- 
tial organization which returns a constant coarse grained 
image of the function cos q, equal in turn to M. Note 
also that a similar coarse grained procedure is adopted 
in [14], with reference to a different case study. In doing 
so, and as a straightforwad consequence of Eq. pTj) . which 
ultimately holds because of the continuum limit hypothe- 
sis, we build up a configuration which returns the sought 
condition: C{x) is a constant function on the circle and 
equal to the magnetization M. The complex spatial or- 
ganization of {p{x),q{x)) yielding to the constant coarse 
grained value of the magnetization M, can be in prin- 
ciple condensed into a local (i.e. dependent on x) one 
particle distribution function f{p,q,x). Clearly, by con- 
struction, any configuration {p(x),q{x)) compatible with 



f{p, q, x), is going to be a solution of the fractional Eq. (130 
However if we aim at building a configuration that shares 
some similarities with the aforementioned QSS, we still 
need stationarity, in the infinite N limit. This property is 
guaranteed if we refer to the stationary distributions that 
were calculated for the HMF in [B]. The recipe goes as 
follows: we pick up q and p values as originating from this 
stationary distribution, so as to ensure that the time evo- 
lution will only consist of a local reshuffling of the actual 
phase space coordinates, without affecting their associated 
coarse grained value. In doing so, we hence obtain a family 
of stationary solutions of the a-HMF model. H 

The microscopic dynamics of the particles in such states 
can be mapped into a pendulum motion and is hence in- 
tegrable. Interestingly, the spatial organization is locally 
scale free. The functions q{x,t) and p{x,t) arc thus "very 
complicated" along the spatial direction, while displaying 
a regular time evolution and no chaos. In order to validate 
this result for a finite size sample on the original lattice we 
proceed as follows. We consider a stationary state of the 
HMF consisting of just one torus with associated mag- 
netization M = 0.5. In pendulum action-angle variable 
we refer hence to a one particle distribution of the type 
/(/, 0) = 5{I — /o)/27r (see for instance [T7] for details re- 
garding the transformation from (p, q) to (/, 9) and vice 
versa). We are focusing our attention on an "individual 
component" belonging to the extended set of a linearly in- 
dependent elements, which define the QSS basis. In order 
to be as close as possible to the stationary state of the a- 
HMF model, we simply distribute randomly on the lattice 
the values picked from such, analytically accessible, dis- 
tribution. The analysis for different values of the number 
of sites and a ~ 0.25, is then performed by monitoring 
the values of the action and shows as expected a trend to- 
wards the stabilization of such a state as N increases (see 
Fig. [T]). We also checked that as expected increasing the 
value of a, which weakens the coupling strength, implies 
inducing a more pronounced destabilization of the state, 
which can be effectively opposed by increasing the num- 
ber of simulated rotors. The solutions here constructed 
are hence stable versus the a-HMF dynamics, provided 
the continuum limit is being performed and so represent a 
consistent analytical prediction for the existence of quasi- 
stationary states, beyond the original HMF setting. 

To further scrutinize the dynamics of the a-HMF model 



^Notc that X corresponds to a point on the circle for (p{x), q{x)), 
while, for for /(p, q, x), it labels the position x of an interval of width 
dx. 

^Wc note that several distinct stationary distributions of the orig- 
inal HMF model are in principle allowed, yielding to identical values 
of AI. It should be hence possible to build up more complex distri- 
butions than those here analyzed, via a combination of such distinct 
profiles assigned to different spatial locations. A stationary distri- 
bution of the a-HMF model can thus correspond to a stationary 
distribution of the HMF model distributed scale free on the circle. 
Alternatively, it can be decomposed in a sum of different stationary 
distributions, which all have the same magnetization and are placed 
in different regions of the circle. Notice that certain restrictions due 
to Eq. (O, may be necessarily imposed at the frontiers. 
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Fig. 1: Top: QSS for the a-HMF model built using a station- 
ary solution of the HMF model corresponding to a one particle 
PDF of the type /(/,6I) = S{I - Io)/2tv (one torus M = 0.5). 
The system is evolved via the a-HMF dynamics with a — 0.25. 
The figure displays the resulting distribution of actions / at 
t = 200, for different sizes. The integration uses the fifth order 
optimal symplectic integrator [15] and a time step 5t = 0.05. 
Bottom: A simple superposition of snapshots relative to dif- 
ferent values of a and at final time t — 200, with A'^ — 2^^. 
a = 0.75, Q = 0.5, a = 0.25, a = 0, correspond respectively 
to the (online) colors red, blue green and black. Clearly the 
thickness of the rings get reduced as a ^ 0, i.e. when the 
a-HMF model tends towards its corresponding HMF limit. 

we turn to direct simulations, starting from out of equi- 
librium initial condition. Our declared goal is to follow 
the system evolution through the violent relaxation pro- 
cess and eventually identify the presence of spontaneoulsy 
emerging QSS for the a-HMF model. In particular, we 
are interested in their microscopic characteristic to make 
a bridge with the analysis developed in the first part of the 
paper. For this purpose we initialize the system in g = 
and assume a Gaussian distribution for the conjugate mo- 
menta p [S] . The system state is monitored by estimating 
the average magnetization amount as a function of the en- 
ergy, see Fig. [2j For energies larger than 0.75 one would 
expect the homogeneous solution to prevail, as dictated 
by the statistical mechanics calculation. However, the sys- 
tem gets confined into a inhomogeneous state, the residual 
time averaged magnetization being large and persistent in 
time. It is therefore tempting to interpret those states as 



Fig. 2: Magnetization vs. energy per particle E/N for different 
a, for a system made of N = 4096 particles. The magnetiza- 
tion values here reported follow from a time average over the 
window 400 < T < 800. 

QSS, and so analyze their associated dynamical features 
in light of the above conclusions. In particular, we expect 
the microscopic dynamics to resemble that of a pendulum, 
bearing some degree of intrinsic regularity. To unravel the 
phase space characteristics we compute the Poincare sec- 
tions, following the recipe in [71[TB] and so visualizing the 
single particle stroboscopic dynamics, with a rate of acqui- 
sition imposed by the self-consistent mean field evolution. 
The averages of C'i and Si refer to the two components 
of the magnetization per site. The Poincare sections are 
drawn by recording the positions pi and qi — ipi in phase 
space each time the equality Mi = M is verified. Re- 
sults for a specific initial conditions are depicted in Fig. [3] 
where one hundred trajectories are retained. The phase 
portrait shares many similarities with that obtained for a 
simple one and a half degree of freedom Hamiltonian (see 
for instance [17)). with many resonances and invariant tori. 
Clearly, and in agreement with the above scenario, a large 
number of particles exhibit regular dynamics. However as 
the nature of phase space reveals, these QSSs are steady 
states of the discrete dynamics, not stationary solutions. 
Nevertheless, we set to analyze the spatial organization 
of the identified steady state to test whether stationary 
state features are present in this configuration. To this 
end, we computed the values of the local magnetization 
M{x, t) = C(x, tY + S(x, tY versus time and estimated 
an individual action, stemming from a Hamiltonian pen- 
dulum, which would give rise to an equation of motion 
formally identical to the Eq.ijSl). Results of the analysis 
are depicted in Fig. |4l One clearly sees that the function 
M{x^t) is homogeneous and presents a dependence on 
thus suggesting that the distribution of q[x^ t) is a solution 
of Eq.([21). The plot of action versus time as depicted in 
Fig. HJ clearly indicates a degree of enhanced spatial com- 
plexity, nearby particles not belonging to the same tori. 
We find in these simulations and in this {N- finite) steady 
state the same distinctive features of the stationary solu- 
tions as depicted earlier. 

To conclude, in this Letter we have shown that station- 
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Fig. 3: Top:"Poincare section" of a QSS for E/N = 1.2, N = 
65536. Initial conditions are Gaussian in p and g = 0. The 
section is computed for 150 < t < 1000 and a — 0.25. Bottom: 
Magnetization versus time. 
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Fig. 4: Top: M (x) = \JC{xY -|- S{xY as a function of time 
,see also Fig. [3] to monitor to appreciate the amplitude of the 
fluctuations as recorded along the time evolution. Bottom: in- 
dividual particle's action as a function of time. The system is 
initialized to get a QSS with ElN = 1.2 and a = 0.25 (see 
Fig. El). Total number of particle is = 8192, x = ijN, with 
1 < i < A''. The actions are those of a pendulum corresponding 
to Eq.Q. M(x) is almost uniform in space, while the individ- 
ual spatial organization is complex. The time evolution of the 
actions appears quite regular. 



ary solutions for the mean field HMF model are as well 
stationary ones of the a-HMF model. Microscopic dy- 
namics in these stationary states is regular, and explicitly 
known. The price to pay for this microscopic regularity in 
time is a complex, self similar, spatial organization corre- 
sponding to the solution of a fractional equation. When 
turning to direct numerical investigations we have identi- 
fied a series of quasi-stationary states, which corresponds 
to steady states. Still, such states share many of the fea- 
tures of their stationary counterparts. The importance 
of these conclusions are manifold: On the one side, we 
confirm that QSS do exist in a generalized non mean field 
setting [TT]. Also, we validate a theoretical method to con- 
struct, from first principle, (out of) equilibrium stationary 
solutions. Finally, the fact that in long range systems sta- 
tionary states (among which we may of course count the 
equilibrium) which display regular microscopic dynamics 
do exist, allows us to dispose of an enormous amount of 
information regarding the intimate dynamics of a system 
frozen in such state. This knowledge can prove crucial 
in bridging the gap between microscopic and macroscopic 
realms in systems with many body interacting elements. 
Interestingly, it could prove useful to investigate the role 
of chaos versus the foundation of statistical mechanics, as 
discussed in pTlfTB] . 
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Appendix. — As mentioned in the core of the letter, 
we shall split the integral (O in L pieces. In this appendix 
we will justify the procedure that eventually leads to eq.( 
llip . To this end and to make the derivation more trans- 
parent we will focus on C(0) and localize the singularity in 
0. The generalization for C\x) is straightforward. Let us 
start by defining the small parameter e— \/L, and write: 
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Consider the intervals labelled with A; > 1, namely those 
where the singularity \\y\\ — is not present. We use the 
regularity of and perform a Taylor expansion of 

the function in the vicinity of j/fc = fc/L to get 
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which implies 
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Assume now that the average J^'j^^^^^ cos q{y)dy is 
constant and equal to M , no matter the value of k. Fur- 
thermore, such a property is assumed to hold as L gets 
larger and larger. Then since cos(y) is bounded by one, 
we can get the following estimate 
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which constraints the Ak term. Exploiting again the fact 
that cos{y) is bound and recalling the definition of C(0) 
we end up with: 
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This relation confirms the adequacy of the approximated 
expression ((TT|) , as it immediately follows by the definition 
of M mentioned above. Furthermore, if we take the limit 
e = 1/L — !> 0, wc get C(0) M, a result which is clearly 
true as long as a < 1. 
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